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Abstract 

Motivated by recent studies of colloidal particles in optical tweezer arrays, we study a two- 
dimensional model of a colloidal suspension in a periodic potential. The particles tend to stay near 
potential minima, approximating a lattice gas. The interparticle interaction, a sum of Yukawa 
terms, features short-range repulsion and attraction at somewhat larger separations, such that two 
particles cannot occupy the same potential well, but occupation of adjacent cells is energetically 
favored. Monte Carlo simulation reveals that the equilibrium system exhibits condensation, as in 
the Ising model/lattice gas with conserved magnetization; the transition appears to be continuous 
at a half occupancy. We study the effect of biased hopping probabilities favoring motion along one 
lattice direction, as might be generated by a steady flow relative to the potential array. This system 
is found to exhibit features of the driven lattice gas: the interface is oriented along the drive, and 
appears to be smooth. A weak drive facilitates ordering of the particles into high- and low-density 
regions, while stronger bias tends to destroy order, and leads to very large energy fluctuations. Our 
results suggest possible realizations of equilibrium and driven lattice gases in a colloidal suspension 
subject to an optical tweezer array. 

PACS numbers: 05.40.-a,82.70.Dd,05.10.Gg,05.50.+q 



Electronic address: potiguar@fisica.ufmg.br 



^Electronic address: dickman@fisica.ufmg.br 



1 



I. INTRODUCTION 



Lattice gas models are an essential tool of statistical mechanics, and have found considerable 
application in the study of phase transitions in adlayers . In the simplest version, 

particles interact only with their nearest neighbors; equilibrium properties of the model are 
well understood The lattice gas is readily extended to the study of time-dependent 
phenomena by defining a particle dynamics (typically via hopping) that obeys detailed bal- 
ance with respect to the Hamiltonian. In driven lattice gas models, also known as driven 
Effusive sterns (CDS), flflQ particles execute biased hopping along one of the principa. 
lattice directions, with periodic boundaries [8|, |9[ , leading to a a nonequilibrium stationary 
state (NESS) with a nonzero current. DDS have been studied extensively as a prototype of 
nonequilibrium phase transitions. Despite their simplicity, DDS exhibit surprising charac- 
teristics jfj for example the increase of the critical temperature with drive strength in the 
case of attractive interactions. A driven fluid in continuous space has also been found to 
possess an unusual phase diagram [lo|. In this case, the drive appears to oppose ordering, 
as reflected in a reduced critical temperature, compared with equilibrium. 
Although much effort has been devoted to studying DDS theoretically and in simulations, 
there are as yet no experimental realizations. Driven colloidal suspensions offer some promise 
in this regard. With the availability of optical-tweezer arrays, the dynamics of colloidal 



particles in an external potential has been investigated intensively [11 
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A 



study of colloidal particles in a periodic potential revealed a variety of dynamic regimes as 



a function of the viscous friction coefficient 



11] . At low friction, the motion at long times 



consists of jumps between adjacent potential minima, resembling that of particles in a lattice 
gas with nearest-neighbor hopping. 

Motivated by this correspondence, we study a two-dimensional model system in continuous 
space, of colloid-like particles in a periodic background potential. The potential is taken 
so that the particles spend most of the time near a potential minimum. The interparticle 
potential, a sum of Yukawa terms, features short-range repulsion and attraction at some- 
what larger separations, such that two particles cannot occupy the same potential well, but 
occupation of adjacent cells is energetically favored. We study the phase behavior of the 
model in Monte Carlo simulations, beginning with the undriven (equilibrium) case, followed 
by an exploration of the effects of driving, and a brief examination of another nonequilibrium 
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version, in which there is no drive, but the external potential is time dependent. 
Although the present model includes, in the interests of computational efficiency, certain 
unrealistic features, we believe that it represents a significant step toward devising a system 
capable of experimental realization, and exhibiting properties characteristic of DDS. Our 
results suggest that it is possible to realize lattice-gas like systems, both equilibrium and 
driven, in a colloidal suspension subject to a periodic external potential. 
In the following section we define the model and simulation method. In Sec. Ill we present 
simulation results, while Sec. IV contains a summary and discussion of open issues, regarding 
both the model system and possible realizations of DDS in experiments on colloids. 



II. MODEL AND SIMULATIONS 



We study a two-dimensional system of N particles interacting via a pairwise potential (r), 
and subject to a periodic external potential V(x, y). The latter is of the form used to model 
an optical tweezer array 15], 



V(x,y) 



Va 



(1) 



where g(x,y) = A [cos(27nr) + cos(27n/) — 2B]. For Vq > the potential minima fall at the 
sites of the integer square lattice (n + 1/2, m + 1/2); parameters A and B together control 
the well depth and the curvature near the potential minimum. The interparticle potential 
is given by the sum of a short-range repulsive term and a slightly longer range attractive 
term, both of Yukawa form : 

-ftir e~ K2r 

u(r) = V 1 V 2 V c , (2) 



for r < r c , and u = for r > r c . Here Kx t 2 are characteristic lengths and V c is a constant 
taken so that u(r c ) = 0. The parameter values used in this study (k\ = 3.30, k 2 = 2.21, 
Vi = 200, V 2 = 90, V Q = 40, A = 5, B = -0.5, and r c = 2.5), lead to a strong short-range 
repulsion that effectively prohibits two particles from occupying the same potential well at 
the temperatures of interest. The minimum of the interparticle potential falls at r = 1, 
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favoring occupation of neighboring cells. (For r = 1, a/2 and 2, one has u = —2.347, -1.35 
and -0.283, respectively.) The ground state energy per particle is -8.396, about 4.2 times 
that of the nearest-neighbor square lattice gas with unit interaction. A crude estimate of the 
critical temperature is then 4.2 times that of the nearest-neighbor lattice gas (T c ~ 0.5673), 
that is, T c ~ 2.4. 

We perform Metropolis Monte Carlo (MC) simulations of the system defined above. The MC 
time step is defined as one attempted move per particle. In each move, a randomly chosen 
particle is subject to a random displacement Ar = (Ax, Ay), with components distributed 
uniformly on the interval [—1,1]. 

As in studies of DDS, the drive takes the form of a force f, such that the work done on a 
particle when it suffers a displacement Ar is f • Ar. As noted, the system is periodic in 
the direction of the drive (in practice, the x direction), so that f cannot be written as the 
gradient of a potential. The acceptance probability for a particle displacement Ar is: 

P(Ar) = min {1, exp [-/3{AE - f • Ar)]} , (3) 



where AE is the change in energy and (3 represents inverse temperature, in units such that 
k B = l- 

The quantities of principal interest are the average energy (E), the order parameter m = 
(m x ,m y ), and the corresponding fluctuations. To define the order parameter we introduce 
lattice-gas variables: <Tij = 1 if the cell centered at (i + 1/2, j + 1/2) is occupied, and is zero 
otherwise. Then an appropriate order parameter for a lattice gas with conserved particle 
density can be defined, for a system of L x L sites, via 

2 

(4) 

with m y given by exchanging indices i and j in the sums. In a disordered phase (m x ) = 
(fflj) = 0. In a half-occupied system, a maximally ordered, isotropic configuration (a square 
of side L/y/2) has m x = m y = 1/2. In a system with periodic boundaries, however, the 
surface energy is smaller in a strip configuration. If the particles occupy all sites the region 
jo < j < jo + L/2, then m x = 1 and m y = 0; the values are exchanged for a strip oriented 
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along the y direction, m = ^m 2 + m 2 characterizes overall ordering. In equilibrium, the 
two orientations are equally likely, but for low temperatures and reasonably large systems, 
the typical time for switching between them is much larger than practical simulation times, 
so it is interesting to characterize the degree of order by the anisotropy parameter Am = 
m> —m < , where m > = max{(m x ), (m y )}, and m < denotes the lesser of the two components. 
Experience with DDS shows that in driven systems m > is always along the drive direction. 

III. NUMERICAL RESULTS 
A. Equilibrium properties 

We study half-occupied systems, N = L 2 /2, in periodic cells of linear size L = 10, 20 and 
30. Our results represent averages over the stationary regime of 10 or more independent 
realizations, each of 10 6 or more time steps. In Fig. [1] we show the mean energy e per 
particle, which exhibits the qualitative behavior typical of lattice gases or fluids with a hard 
core and short-range attraction. Evidence of a phase transition is seen in the progressive 
growth, with system size, of the peak in the specific heat per particle, c = vai(E)/NT 2 , 
shown in Fig. 2. The specific heat maxima fall at temperatures T m = 2.1, 2.5 and 2.63 
for sizes L = 10, 20 and 30, respectively. These data are fit very well by the expression 
T m = T c + Const./ L, yielding the estimate T c = 2.90(1) for the transition temperature. 
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FIG. 1: Energy per particle versus temperature for L = 10, 20 and 30 (upper to lower). 
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FIG. 2: Specific heat per particle versus temperature for L = 10, 20 and 30, in ascending order. 

In Fig. 3 we show the order parameter m and the anisotropy Am for L = 30. These results 
are consistent with a transition near T = 2.9. Of note is the much smaller finite size effect 
in Am as compared with m in the high-temperature regime, as the fluctuations in the two 
components of m are nearly equal. The inset shows that x = -rp-var(m) exhibits a well 
defined maximum in the vicinity of the transition. 

Based on the present results, we can conclude that the half-filled system exhibits an (ap- 
parently continuous) order-disorder transition near a temperature of 2.9. Further details on 
the nature of the transition (which from symmetry considerations should belong to the Ising 
model universality class), must await more extensive studies using larger systems. 
To close this section we show a representative configuration (Fig. 4) for L = 30 at tempera- 
ture T = 2. Here the system has separated into distinct high- and low-density phases. The 
latter is quite dilute, while the former clearly reflects the periodic potential, and possesses a 
low density of vacant sites. The interface between the high- and low-density regions is quite 
rough. 

B. Driven system 

We determined the stationary properties for various drive strengths /, for system sizes 
L = 20 and L = 30, following the procedure described above. The mean energy is plotted 
versus temperature in Fig. for various drive intensities. For the relatively weak drive 
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FIG. 3: Magnetization m (filled symbols) and anisotropy Am (open symbols) versus temperature 
for L = 30. Inset: scaled variance of the order parameter. 




FIG. 4: Typical configuration, L = 30, T = 2. 

(/ = 1), the energy is just slightly higher than in equilibrium. For the larger drives we 
observe a substantial increase in energy. 

In Fig. [6] we compare the scaled energy variance c =var(E)/NT 2 in equilibrium and under 
drive / = 1. Although this quantity does not represent the specific heat for f > 0, it is 
nevertheless reasonable to suppose that a singularity in c (or a sharp peak, in a finite system), 
marks a phase transition. It is therefore interesting to note that the / = 1 data exhibit a 
sharper peak (and at a slightly higher temperature), than in equilibrium. This suggests 
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FIG. 6: c =\ai{E) /NT 2 versus T for drive / = 1 (upper curve) and in equilibrium (lower). Upper 
inset: comparison, for / = 1, of the scaled variance c (points) and de/dT (smooth curve). Lower 
inset: scaled variance for / = 1 and L = 20 (open symbols) and L = 30 (filled symbols). 

that a weak drive facilitates ordering. For / = 1, the scaled variance (upper inset) follows 
the same trends as, but is generally greater than, the specific heat, de/dT. (In this case 
we obtain de/dT via numerical differentiation of a polynomial fit to the energy data.) This 
suggests that the drive causes fluctuations beyond those generated by thermal mechanisms. 
(That is, the drive realizes work on the system but is not included in the energy balance.) 
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The dependence of var(e) on system size (lower inset) is qualitatively similar to that found 
in equilibrium; we estimate T c = 2.85(5) for / = 1. 

The plot of the order parameter (Fig. [7]) confirms that a weak drive (/ = 1, 2) enhances 
ordering, whereas a stronger one inhibits it. These data suggest a transition temperature of 
T c « 2, for / = 7. 
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FIG. 7: Order parameter versus temperature for L = 30 and drive / = (open squares), / = 1 
(+), / = 2 (open circles), / = 4 (diamonds), and / = 7 (filled squares). 

The effect of increasing drive at fixed temperature is shown in Figs. [8] (energy and its 
variance) and [9] (order parameter). The energy increases slowly with / for a weak drive, 
rapidly for intermediate drive strength, and then exhibits a steady, more gradual growth for 
/ > 8 or so. The energy variance shows a marked peak (near / = 4 for T = 2.5, and / = 6 
for T = 2.0), which appears to be associated with destruction of an ordered arrangement. 
The amplitude of this peak is much larger than in equilibrium. For larger drives, var(e) 
increases steadily with /, as the drive forces particles out of the periodic potential minima. 
Fig. [9] again shows that a weak drive enhances ordering. The order parameter reaches a 
maximum near / = 1 - 2, and then decays rapidly when the drive is increased further, and 
the driving force begins to dominate interparticle attraction. 

The stationary current density j, defined as the mean displacement (Az) per site and unit 
time, is plotted as a function of drive in Fig. [TOj For the range of parameters studied here, 
j is an increasing function of both / and T. In the disordered phase (upper set of points in 
Fig. [TOl) . the current grows linearly with / for small /, and then shows signs of saturating 
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FIG. 9: Order parameter components m > (upper curves) and m < (lower curves) versus drive / for 
T = 2.5 (open symbols), and T = 2.0 (filled symbols), L = 30. Inset: m = \J m> + m< versus /. 

at larger values of the drive. At lower temperatures and weak drives, such that the system 
is ordered, the current is severely reduced, but it takes values comparable to those at higher 
temperature once / is large enough to disorder the system. The latter event is signaled by 
a sharp peak in the variance of the energy (Fig. [8]); the singularity (if any) in the current is 
much weaker. 

Our results suggest that the particles are organized into dense stripes along the field direc- 
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FIG. 10: Current density j versus drive / for T = 3.5 (filled symbols), and T = 2.0 (open symbols), 
L = 30. 

tion, as in the driven lattice gas. The configurations shown in Figs. [11] (T = 2, / = 1), 
and fT2l (T = 1.1, / = 1), confirm this expectation. Comparing Figs. HI (/ = 0) and fTTI 
(/ = 0), the chief difference is the smoother interface found in the driven system. Fig. [T2l 
shows a configuration well below the transition temperature, in which two separate stripes 
have emerged. As in the driven lattice gas, this appears to be a long-lived metastable state 
at low temperature [7| . 




FIG. 11: Typical configuration, L = 30, / = 1.0, T = 2.0. The drive is directed to the right. 
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FIG. 12: Typical configuration, L = 30, / = 1.0, T = 1.1. 



C. Moving background potential 

We performed some preliminary studies of another method for perturbing the system out 
of equilibrium. In this case there is no drive, but the external potential is time dependent, 
given by V(x — vt, y), with the function V(x, y) as in Eq. ([T]). [Since the potential is periodic, 
the first argument of V is effectively x — vt (mod 1).] The particles are dragged along by 
the moving potential array, as in optical peristalsis 17]. 

In these studies we use L = 30; the background potential amplitude Vq is reduced from 
40 to 10, and the trial particle displacement Ax is uniform on [—.5, .5] (similarly for Ay). 
Otherwise the parameters are as in the studies reported above. For velocities v < 0.02, 
there is anisotropic, lattice-gas like ordering for T < T c . (The critical temperature is about 
2.5 for v = 0.01.) For larger velocities, v = 0.05 or 0.1, the critical temperature is reduced 
considerably (T c ~ 1.8 for v = 0.1). The interface between phases is again oriented along the 
x direction, but for v > 0.05 the particle configuration (Fig. [T3]) shows no hint of periodic 
structure. 
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FIG. 13: Typical configuration under moving background potential, L = 30, v = 0.1, T = 1.5. 
IV. CONCLUSIONS 

We study a two dimensional model of colloidal particles subject to a periodic external po- 
tential and an interaction that favors ordering in a lattice-gas-like arrangement. Stationary 
properties of the model are determined via Monte Carlo simulation. In equilibrium the 
system appears to undergo a second order phase transition, as reflected in the specific heat 
and the order parameter. When the short-range hopping dynamics is biased by a drive, 
we observe anisotropic phase separation, i.e., the interface is oriented along the drive, and 
becomes smoother, as found in studies of the driven lattice gas [3] . A weak drive facilitates 
ordering, as in the driven lattice gas. A strong drive, by contrast, tends to destroy order, 
and provokes very large energy fluctuations. Preliminary analysis indicates that a colloidal 
system under a moving background potential may offer another realization of anisotropic 
phase separation. 

Our results suggest that a system of particles interacting via a potential featuring short- 
range repulsion and slightly longer range attraction, in the presence of a periodic external 
potential with a lattice constant matched to the interparticle potential minimum, provides 
a good candidate for realization of both equilibrium and driven lattice gases. Although the 
physics of such a continuous space system is richer than the corresponding lattice model, 
there is good reason to expect that similar scaling properties will emerge, in a suitable range 
of parameter values. 
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While the equilibrium model could in principle be realized using a planar optical tweezer 
array to provide the external potential, the driven system must, by nature, have periodic 
boundary conditions along the drive direction. This suggests a cylindrical geometry, with the 
driving force provided by Stokes drag on the particles in steady Couette flow. (This scheme, 
without the optical tweezer array, could be used to realize the type of system studied in 
Ref. [10].) Alternatively, one could impose relative rotational motion between the colloidal 
suspension and the tweezer array. We defer analysis of the feasibility of such a setup to 
future work. 

Many additional aspects of the system remain to be explored theoretically and in simulations. 
While the present study uses (because of its greater efficiency) Monte Carlo simulation to 
map out equilibrium properties and nonequilibrium steady states, the driven system should 
be investigated via Langevin dynamics. Aside from affording a better account of nonequi- 
librium stationary and time-dependent properties, this method would allow for inclusion 
of hydrodynamic interactions between the particles and with container walls. Preliminary 
studies of the model system studied here using the Langevin equation in fact yield qualita- 
tively similar results to those reported above [l9|]. It would also be interesting to remove 
the periodic boundaries in the direction perpendicular to the drive, and the restriction to 
two dimensions. More detailed characterization of the phase transitions exhibited by this 
system, both in and out of equilibrium, using large-scale simulations, are planned for future 
work. 
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